/* -------------------------------------------------------------------------- */
/* UMFPACK Copyright (c) 2005-2012 by Timothy A. Davis,
/* Modified by Min Li, Sep, 2013
/* http://www.suitesparse.com. All Rights Reserved.  See ../Doc/License       */
/* for License.                                                               */
/* -------------------------------------------------------------------------- */

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "umfpack.h"

//int    n = 5 ;
//int    Ap [ ] = {0, 2, 5, 9, 10, 12} ;
//int    Ai [ ] = { 0,  1,  0,   2,  4,  1,  2,  3,   4,  2,  1,  4} ;
//double Ax [ ] = {2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1.} ;
//double b [ ] = {8., 45., -3., 3., 19.} ;
//double x [5] ;

#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#ifndef TRUE
#define TRUE (1)
#endif
#ifndef FALSE
#define FALSE (0)
#endif


static double resid
(
    int transpose,
    int Ap [ ],
    int Ai [ ],
    double Ax [ ],
    int n,
    double *b,
    double *x
)
{
    int i, j, p ;
    double norm ;
    double r[n];
    
    for (i = 0 ; i < n ; i++)
    {
            r [i] = -b [i] ;
    }
    if (transpose)
    {
	for (j = 0 ; j < n ; j++)
	{
	    for (p = Ap [j] ; p < Ap [j+1] ; p++)
	    {
		i = Ai [p] ;
		r [j] += Ax [p] * x [i] ;
	    }
	}
    }
    else
    {
	for (j = 0 ; j < n ; j++)
	{
	    for (p = Ap [j] ; p < Ap [j+1] ; p++)
	    {
		i = Ai [p] ;
		r [i] += Ax [p] * x [j] ;
	    }
	}
    }
    norm = 0. ;
    for (i = 0 ; i < n ; i++)
    {
	norm = MAX (ABS (r [i]), norm) ;
    }
    return (norm) ;
}


int main (int argc, char* argv[])
{
        char line[65536];
        int i, n, nz, nz1, status;
        int *A_row, *A_col, *Ap, *Ai;
        double *A_val, *Ax, *B, *X;
        void *Symbolic, *Numeric ;
        
        if(argc != 6){
                printf("Usage: ./umfpack_HSTA [A_info] [A_row] [A_col] [A_val] [B]\n");
                exit(1);
        }
                
        double t[2], rnorm;
        double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
        
        /* ---------------------------------------------------------------------- */
        /* initializations */
        /* ---------------------------------------------------------------------- */
        
        umfpack_tic (t) ;
        
        printf ("\nUMFPACK V%d.%d (%s) demo: _di_ version\n",
                UMFPACK_MAIN_VERSION, UMFPACK_SUB_VERSION, UMFPACK_DATE) ;
        
        /* get the default control parameters */
        umfpack_di_defaults (Control) ;
        
        /* change the default print level for this demo */
        /* (otherwise, nothing will print) */
        Control [UMFPACK_PRL] = 6 ;
        
        /* print the control parameters */
        umfpack_di_report_control (Control) ;
    
        
        /* ---------------------------------------------------------------------- */
        /* read in the input files */
        /* ---------------------------------------------------------------------- */
        FILE* f_A_info;
        f_A_info = fopen(argv[1], "r");
        assert(f_A_info != NULL);
        fgets(line, 65536, f_A_info);
        n = atoi(line);
        printf("i: Dimention %d\n", n);
        
        fgets(line, 65536, f_A_info);
        nz = atoi(line);
        printf("i: Non-zero %d\n", nz);
        fclose(f_A_info);
        
        //malloc arrays 
        nz1 = (nz > 1) ? nz : 1;	/* ensure arrays are not of size zero. */
        A_row = (int *) malloc ((nz) * sizeof (int)) ;
        A_col = (int *) malloc ((nz) * sizeof (int)) ;
        A_val = (double *) malloc ((nz) * sizeof (double)) ;
        B = (double *) malloc (n * sizeof (double)) ;
        Ap = (int *) malloc ((n+1) * sizeof (int)) ;
        Ai = (int *) malloc (nz1 * sizeof (int)) ;
        Ax = (double *) malloc (nz1 * sizeof (double)) ;
        X = (double *) malloc (n * sizeof (double)) ;
        
        if (!A_row || !A_col || !Ap || !Ai || !Ax){
                error ("out of memory") ;
        }

        FILE* f_A_row;
        f_A_row = fopen(argv[2], "r");
        assert(f_A_row != NULL);
        for(i = 0; i < nz1; i++){
                fgets(line, 65536, f_A_row);
                A_row[i] = atoi(line);
                //printf("r: %d\n", A_row[i]);
        }
        
        FILE* f_A_col;
        f_A_col = fopen(argv[3], "r");
        assert(f_A_col != NULL);
        for(i = 0; i < nz1; i++){
                fgets(line, 65536, f_A_col);
                A_col[i] = atoi(line);
                //printf("c: %d\n", A_col[i]);
        }
        
        FILE* f_A_val;
        f_A_val = fopen(argv[4], "r");
        assert(f_A_val != NULL);
        for(i = 0; i < nz1; i++){
                fgets(line, 65536, f_A_val);
                A_val[i] = atof(line);
                //printf("v: %.5f\n", A_val[i]);
        }
        
        FILE* f_B;
        f_B = fopen(argv[5], "r");
        assert(f_B != NULL);
        for(i = 0; i < n; i++){
                fgets(line, 65536, f_B);
                B[i] = atof(line);
                //printf("b: %.5f\n", B[i]);
        }
        
        /* print the right-hand-side */
        printf ("\nb: ") ;
        (void) umfpack_di_report_vector (n, B, Control) ;
    
        /* print the triplet form of the matrix */
        printf ("\nA: ") ;
        (void) umfpack_di_report_triplet (n, n, nz, A_row, A_col, A_val, Control) ;
        
        status = umfpack_di_triplet_to_col (n, n, nz, A_row, A_col, A_val,
                                            Ap, Ai, Ax, (int *) NULL) ;
        if (status < 0){
                //umfpack_di_report_status (Control, status) ;
                error ("umfpack_di_triplet_to_col failed") ;
        }
        
        /* print the column-form of A */
        printf ("\nA: ") ;
        //(void) umfpack_di_report_matrix (n, n, Ap, Ai, Ax, 1, Control) ;
        
        /* ---------------------------------------------------------------------- */
        /* symbolic factorization */
        /* ---------------------------------------------------------------------- */
        double *null = (double *) NULL ;
        status = umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, Control, Info) ;
        if (status < 0){
                printf("Symbolic factorization fail\n");
                error ("umfpack_di_symbolic failed") ;
        }
        
        /* print the symbolic factorization */
        printf ("\nSymbolic factorization of A: ") ;
        //status = umfpack_di_report_symbolic (Symbolic, Control) ;
        if (status < 0){
                error ("umfpack_di_symbolic report failed") ;
        }
        
        /* print the numeric factorization */
        printf ("\nNumeric factorization of A: ") ;
        status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, Control, Info) ;
        if (status < 0){
                umfpack_di_report_info (Control, Info) ;
                umfpack_di_report_status (Control, status) ;
                printf("Numeric factorization fail\n");
                error ("umfpack_di_numeric failed") ;
        }
        
        //status = umfpack_di_report_numeric (Numeric, Control) ;
        if (status < 0){
                error ("umfpack_di_numeric report failed") ;
        }
        
        umfpack_di_free_symbolic (&Symbolic) ;
        status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, X, B, Numeric, Control, Info) ;
        umfpack_di_report_info (Control, Info) ;
        umfpack_di_report_status (Control, status) ;
        if (status < 0){
                printf("Solving fail\n");
                error ("umfpack_di_solve failed") ;
        }
        
        umfpack_di_free_numeric (&Numeric) ;
        
        //print results
        printf ("\nx (solution of Ax=b): ") ;
        printf("Printing results...\n");
        status = umfpack_di_report_vector (n, X, Control) ;
        if (status < 0){
                error ("umfpack_di_report_vec failed") ;
        }
        rnorm = resid (FALSE, Ap, Ai, Ax, n, B, X) ;
        printf ("norm is %g\n", rnorm) ;

        char *str = (char*)malloc((strlen(argv[1])+1) * sizeof(char));
        strcpy(str, argv[1]);
        
        for(i = strlen(str) - 1; i >= 0; i--)
                if(str[i] == '.')
                        break;
        char filename[60];
        strncpy(filename, str, i);
        filename[i] = '\0';
        strcat(filename, ".out");
        printf("Output file name is %s\n", filename);
        
        FILE* f_sol;
        f_sol = fopen(filename, "w+");
        assert(f_sol != NULL);
        for (i = 0 ; i < n ; i++) fprintf (f_sol, "%g\n", X[i]) ;
        fclose(f_sol);
        
        return (0) ;
}

